8.1 KB · ⇊ 3 parts · ⤓ Micron 8.1 KB · ⤓ HTML 56.4 KB

top

Das implizite Trapez-Verfahren ist ein Verfahren zur numerischen Lösung eines Anfangswertproblems

y ′ ( t ) = f ( t , y ( t ) ) , y ( t 0 ) = y 0 {\displaystyle y'(t)=f\left(t,y(t)\right),\quad y(t_{0})=y_{0}}

Es lässt sich sowohl den Runge-Kutta-Verfahren als auch den Adams-Moulton-Verfahren zuordnen. Das Trapezverfahren ist A-stabil mit der Besonderheit, dass für die Schwingungsgleichung y ′ = i α α y {\displaystyle y'=\mathrm {i} \alpha y} kein Amplitudenfehler auftrittcite-ref-1[1]. Das Verfahren lässt sich aus der Trapezregel herleiten:

y n + 1 = y n + h 2 ( f n + 1 + f n ) {\displaystyle y_{n+1}=y_{n}+{\frac {h}{2}}(f_{n+1}+f_{n})}

mit

f n := f ( t n , y n ) . {\displaystyle f_{n}\ :=f(t_{n},y_{n}).}

• Herleitung

• Lösungsmethode

• Stabilität

• Schrittweite h

• Literatur

• Einzelnachweise

Für die Herleitung von Einschrittverfahren wird das Anfangswertproblem meist in der zu ihr äquivalenten Integralgleichung umgeformtcite-ref-2[2]

y ˙ ˙ = f ( t , y ) , y ( t 0 ) = y 0 ⟺ ⟺ y ( t ) = y 0 + ∫ ∫ t 0 t f ( s , y ( s ) ) d s . {\displaystyle {\begin{aligned}{\dot {y}}&=f(t,y)\,,\qquad y(t_{0})=y_{0}\\\Longleftrightarrow \quad y(t)&=y_{0}+\int _{t_{0}}^{t}f(s,y(s))\,\mathrm {d} s\,.\end{aligned}}}

Nun besteht die Idee bei der impliziten Trapez-Methode eine simple Quadraturformel für das Integral zu benutzen: die Trapezregel. Man approximiert in jedem k {\displaystyle k} -ten Schritt den Integranden wie folgt

∫ ∫ t k t k + 1 f ( s , y ( s ) ) d s ≈ ≈ h 2 ( f ( t k , y k ) + f ( t k + 1 , y k + 1 ) ) . {\displaystyle \int _{t_{k}}^{t_{k+1}}f(s,y(s))\,\mathrm {d} s\approx {\frac {h}{2}}{\Big (}f(t_{k},y_{k})+f(t_{k+1},y_{k+1}){\Big )}\,.}

Zusammen ergibt dies die Trapez-Methodecite-ref-3[3]

y ( t k + 1 ) = y ( t k ) + ∫ ∫ t k t k + 1 f ( s , y ( s ) ) d s ≈ ≈ y k + h 2 ( f ( t k , y k ) + f ( t k + 1 , y k + 1 ) ) =: y k + 1 . {\displaystyle y(t_{k+1})=y(t_{k})+\int _{t_{k}}^{t_{k+1}}f(s,y(s))\,\mathrm {d} s\approx y_{k}+{\frac {h}{2}}{\Big (}f(t_{k},y_{k})+f(t_{k+1},y_{k+1}){\Big )}=:y_{k+1}\,.}

Zur Lösung dieses, in der Regel nichtlinearen, Gleichungssystems können verschiedene numerische Verfahren genutzt werden. Für das quadratisch konvergente Newton-Verfahren ergibt sich konkret:

y n + 1 ( k + 1 ) = y n + 1 ( k ) − − ( I − − h 2 ∂ ∂ f n + 1 ( k ) ∂ ∂ y n + 1 ( k ) ) − − 1 ( y n + 1 ( k ) − − y n − − h 2 ( f n + 1 ( k ) + f n ) ) . {\displaystyle y_{n+1}^{(k+1)}=y_{n+1}^{(k)}-\left(I-{\frac {h}{2}}{\frac {\partial f_{n+1}^{(k)}}{\partial y_{n+1}^{(k)}}}\right)^{-1}\left(y_{n+1}^{(k)}-y_{n}-{\frac {h}{2}}(f_{n+1}^{(k)}+f_{n})\right).}

Man erhält also ein lineares Gleichungssystem

( I − − h 2 J ( k ) ) y n + 1 ( k + 1 ) = − − h 2 J ( k ) y n + 1 ( k ) + y n + h 2 ( f n + 1 ( k ) + f n ) , {\displaystyle (I-{\frac {h}{2}}J^{(k)})y_{n+1}^{(k+1)}=-{\frac {h}{2}}J^{(k)}y_{n+1}^{(k)}+y_{n}+{\frac {h}{2}}(f_{n+1}^{(k)}+f_{n}),}

wobei J die Jacobi-Matrix

J ( k ) := ( ∂ ∂ f ∂ ∂ y ) n + 1 ( k ) {\displaystyle J^{(k)}:=\left({\frac {\partial f}{\partial y}}\right)_{n+1}^{(k)}} ,

I {\displaystyle I} die Einheitsmatrix und k {\displaystyle k} der Iterationsschritt ist.

Mit der Testgleichung y ′ ( t ) = λ λ y ( t ) {\displaystyle y'(t)=\lambda y(t)} bekommt man die Stabilitätsfunktion

R ( z ) = 2 + z 2 − − z , z = h λ λ ∈ ∈ C . {\displaystyle R(z)={\frac {2+z}{2-z}},\quad z=h\lambda \in \mathbb {C} .}

Auf der imaginären Achse z = i η η {\displaystyle z=i\eta } gilt | R ( i η η ) | = 1 {\displaystyle |R(i\eta )|=1} , daher ist die Trapezmethode A-stabil.

Die (variable) Schrittweite kann aus folgender Beziehung berechnet werden:

| R ( h λ λ ) e h λ λ − − 1 | = δ δ {\displaystyle \left\vert {\frac {R(h\lambda )}{\mathrm {e} ^{h\lambda }}}-1\right\vert =\delta } ;

δ δ {\displaystyle \delta } bezeichnet den zugelassenen lokalen Diskretisierungsfehler. Der Ansatz y n + 1 = y n + h 2 ( f n + 1 + f n ) =: R ( h λ λ ) y n {\displaystyle y_{n+1}=y_{n}+{\frac {h}{2}}(f_{n+1}+f_{n})=:R(h\lambda )y_{n}} liefert für die implizite Trapez-Methode

R ( h λ λ ) = 2 + h λ λ 2 − − h λ λ {\displaystyle R(h\lambda )={\frac {2+h\lambda }{2-h\lambda }}} .

Dabei ist λ λ := max j | λ λ j | {\displaystyle \lambda \,:=\max _{j}{|\lambda _{j}|}} der Betrag des betragsmäßig größten Eigenwerts der Jacobi-Matrix (Spektralradius). Die numerische Bestimmung der Eigenwerte ist sehr zeitaufwendig; für den Zweck der Schrittweitenberechnung ist es im Allgemeinen ausreichend die Gesamtnorm λ λ = N ⋅ ⋅ max i , j | a i j | {\displaystyle \lambda =N\cdot \max _{i,j}|a_{ij}|} heranzuziehen, die immer größer oder gleich der Spektralnorm ist. N ist der Rang der Jacobi-Matrix und a i j {\displaystyle a_{ij}} deren Elemente.

• Hans R. Schwarz, Norbert Köckler: Numerische Mathematik. 5. Auflage, Teubner, Stuttgart 2004, ISBN 3-519-42960-8, S. 343.

cite-note-11. ↑ M. Kloker: Numerische Löser (Zeitintegrationsverfahren) für die Gewöhnliche Modelldifferentialgleichung y'=αy (PDF; 2,2 MB), Universität Stuttgart, 1996

cite-note-22. ↑ Reusken, Arnold.: Numerik für Ingenieure und Naturwissenschaftler. Springer, Berlin 2006, ISBN 3-540-25544-3, S. 378.

cite-note-33. ↑ Reusken, Arnold.: Numerik für Ingenieure und Naturwissenschaftler. Springer, Berlin 2006, ISBN 3-540-25544-3, S. 383.